library(ggordiplots)
library(microViz)
library("vegan")
library(phyloseqCompanion)
library(ggplot2)
library(Rmisc)
library(ggsignif)
library(ggpubr)
library(dplyr)
library("gridExtra")
library(microbiomeMarker)
library(tidyr)
library(FUNGuildR)
library(betareg)
library(SciViews)
library(MuMIn)
library(MASS)
library(mblm)
library(lme4)
library(splines)
library(lmerTest)

# set scientific notation off
options(scipen=999)

setwd("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code")

# load data
load("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code/RA otu by guild and exploration type.RData")

# load phyloseq object and exploration type and guild data
load("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code/RA otu by guild and exploration type_new.RData")

# load newest data
load("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing/Code/functional type by sample.RData")

###########################################################################################################################
# test differences in the relative abundance of exploration types by Hansen close bank
###########################################################################################################################

# remove right bank distances near decomposing salmon carcass?
exploration.guild.RA <- exploration.guild.RA[-c(42,70,81),]

# or don't include S4 SR 1
exploration.guild.RA <- exploration.guild.RA[-c(70,81),]

# close bank (1-6m)
hansen.left.close <- exploration.guild.RA$short.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SL")]
hansen.right.close <- exploration.guild.RA$short.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SR")]
t.test(hansen.left.close, hansen.right.close)
one.way <- aov(short.distance ~ Side, data=exploration.guild.RA[which(exploration.guild.RA$Stream== "Hansen" & exploration.guild.RA$Distance < 7),])
summary(one.way)
wilcox.test(hansen.left.close, hansen.right.close, exact=F)

hansen.left.close <- exploration.guild.RA$medium.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SL")]
hansen.right.close <- exploration.guild.RA$medium.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SR")]
t.test(hansen.left.close, hansen.right.close)
wilcox.test(hansen.left.close, hansen.right.close, exact=F)
one.way <- aov(medium.distance ~ Side, data=exploration.guild.RA[which(exploration.guild.RA$Stream== "Hansen" & exploration.guild.RA$Distance < 7),])

# medium distance higher on right bank when removing decomposing carcasses S4 SR 1-6

hansen.left.close <- exploration.guild.RA$long.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SL")]
hansen.right.close <- exploration.guild.RA$long.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SR")]
t.test(hansen.left.close, hansen.right.close)
wilcox.test(hansen.left.close, hansen.right.close, exact=F)

# close bank (3-6m)
hansen.left.close <- exploration.guild.RA$short.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SL")]
hansen.right.close <- exploration.guild.RA$short.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SR")]
t.test(hansen.left.close, hansen.right.close)
wilcox.test(hansen.left.close, hansen.right.close, exact=F)

hansen.left.close <- exploration.guild.RA$medium.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SL")]
hansen.right.close <- exploration.guild.RA$medium.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SR")]
t.test(hansen.left.close, hansen.right.close)
wilcox.test(hansen.left.close, hansen.right.close, exact=F)

hansen.left.close <- exploration.guild.RA$long.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SL")]
hansen.right.close <- exploration.guild.RA$long.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SR")]
t.test(hansen.left.close, hansen.right.close)
wilcox.test(hansen.left.close, hansen.right.close, exact=F)

# no difference except MAYBE long distance between 3-6 m banks, but only for full data (not for when remove S4 SR 1-6)

# test by full Hansen bank
hansen.left <- exploration.guild.RA$short.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SL")]
hansen.right <- exploration.guild.RA$short.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SR")]
t.test(hansen.left, hansen.right)
wilcox.test(hansen.left, hansen.right, exact=F)

hansen.left.close <- exploration.guild.RA$medium.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SL")]
hansen.right.close <- exploration.guild.RA$medium.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SR")]
t.test(hansen.left, hansen.right)
wilcox.test(hansen.left, hansen.right, exact=F)

hansen.left.close <- exploration.guild.RA$long.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SL")]
hansen.right.close <- exploration.guild.RA$long.distance[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SR")]
t.test(hansen.left, hansen.right)
wilcox.test(hansen.left, hansen.right, exact=F)

# no difference in exploration types between Hansen full bank

#########################################
# test differences in the relative abundance of fungal guilds by Hansen close bank
#########################################

# close bank (1-6m)
left <- exploration.guild.RA$saprotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SL")]
right <- exploration.guild.RA$saprotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SR")]
t.test(left, right)
wilcox.test(left, right, exact=F)

left <- exploration.guild.RA$symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SL")]
right <- exploration.guild.RA$symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SR")]
t.test(left, right)
wilcox.test(left, right, exact=F)

# maybe symbiotrophs different

left <- exploration.guild.RA$saprotroph.symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SL")]
right <- exploration.guild.RA$saprotroph.symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Side=="SR")]
t.test(left, right)
wilcox.test(left, right, exact=F)

# maybe saprotroph.symbiotroph different

# close bank (3-6m)
left <- exploration.guild.RA$saprotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SL")]
right <- exploration.guild.RA$saprotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SR")]
t.test(left, right)
wilcox.test(left, right, exact=F)

left <- exploration.guild.RA$symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SL")]
right <- exploration.guild.RA$symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SR")]
t.test(left, right)
wilcox.test(left, right, exact=F)

left <- exploration.guild.RA$saprotroph.symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SL")]
right <- exploration.guild.RA$saprotroph.symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1 & exploration.guild.RA$Side=="SR")]
t.test(left, right)
wilcox.test(left, right, exact=F)

# full bank
left <- exploration.guild.RA$saprotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SL")]
right <- exploration.guild.RA$saprotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SR")]
t.test(left, right)
wilcox.test(left, right, exact=F)

left <- exploration.guild.RA$symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SL")]
right <- exploration.guild.RA$symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SR")]
t.test(left, right)
wilcox.test(left, right, exact=F)

left <- exploration.guild.RA$saprotroph.symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SL")]
right <- exploration.guild.RA$saprotroph.symbiotroph[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SR")]
t.test(left, right)
wilcox.test(left, right, exact=F)

# plot

# choose whether full bank or close bank at Hansen
hansen.types <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen"),]
hansen.types.close <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7),]
hansen.types.real.close <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Distance < 7 & exploration.guild.RA$Distance > 1),]

# plot long distance from 3-6 m 
tgc <- summarySE(hansen.types.real.close, measurevar="long.distance", groupvars=c("Side"))
tgc

long.distance <- ggplot(tgc, aes(Side, long.distance, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Long-distance types")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=long.distance-se, ymax=long.distance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14),
        legend.position="none")

# plot symbiotrophs from 1-6 m
tgc <- summarySE(hansen.types.close, measurevar="symbiotroph", groupvars=c("Side"))
tgc

symbio <- ggplot(tgc, aes(Side, symbiotroph, fill=Side)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Symbiotrophs")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=symbiotroph-se, ymax=symbiotroph+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14),
        legend.position="none")

p <- grid.arrange(long.distance, symbio, nrow=1, ncol=2)
p
ggsave(plot=p, width=10, height=6, dpi=300, filename = "Fig2.jpg")

###########################################################################################################################
# test differences in relative abundance of fungal guilds by C:N at all streams
###########################################################################################################################

# for beta regression, need to remove any zeros and replace them with almost zeros since beta regression only takes (0,1) in dependent variable
exploration.guild.RA$saprotroph_trans <- c(1:95)
exploration.guild.RA$symbiotroph_trans <- c(1:95)
exploration.guild.RA$medium.distance_trans <- c(1:95)
exploration.guild.RA$long.distance_trans <- c(1:95)

for (i in 1:95){
  if (exploration.guild.RA$saprotroph[i] == 0.0000000000){
    exploration.guild.RA$saprotroph_trans[i] <- 0.000000000000000000001
  }
  else {exploration.guild.RA$saprotroph_trans[i] <- exploration.guild.RA$saprotroph[i]}
}

for (i in 1:95){
  if (exploration.guild.RA$symbiotroph[i] == 0.0000000000){
    exploration.guild.RA$symbiotroph_trans[i] <- 0.000000000000000000001
  }
  else {exploration.guild.RA$symbiotroph_trans[i] <- exploration.guild.RA$symbiotroph[i]}
}

for (i in 1:95){
  if (exploration.guild.RA$medium.distance[i] == 0.0000000000){
    exploration.guild.RA$medium.distance_trans[i] <- 0.000000000000000000001
  }
  else {exploration.guild.RA$medium.distance_trans[i] <- exploration.guild.RA$medium.distance[i]}
}

for (i in 1:95){
  if (exploration.guild.RA$long.distance[i] == 0.0000000000){
    exploration.guild.RA$long.distance_trans[i] <- 0.000000000000000000001
  }
  else {exploration.guild.RA$long.distance_trans[i] <- exploration.guild.RA$long.distance[i]}
}

#########################################
### FULL MODEL BY GUILD ###
#########################################

# remove right bank distances near decomposing salmon carcass?
exploration.guild.RA <- exploration.guild.RA[-c(42,70,81),]

# or keep S4 SR 1 (but remove S4 SR 3 and S4 SR 6)
exploration.guild.RA <- exploration.guild.RA[-c(70,81),]

# hansen only data
hansen <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Hansen"),]

# examine data by stream
ggplot(exploration.guild.RA,aes(x=long.diversity, y=Stream)) + geom_boxplot(alpha=0.2) + stat_summary(fun=mean, geom="point")

# after examining relative abundance of guilds by stream, concluded that stream doesn't have to be included as a random effect

# non parametric kendall-theill regresson
mod1 <- mblm(saprotroph ~ C.N, data=exploration.guild.RA)
mod2 <- mblm(symbiotroph ~ C.N, data=exploration.guild.RA)
mod3 <- mblm(saprotroph.symbiotroph ~ C.N, data=exploration.guild.RA)
summary(mod1)
summary(mod2)
summary(mod3)

# saprotrophs decreased with C:N, symbiotrophs increased with C:N, and NS for combined guild (p-0.06, +)

# LM and GLM
mod1 <- lm(saprotroph ~ C.N, data=exploration.guild.RA)
mod2 <- lm(symbiotroph ~ C.N, data=exploration.guild.RA)
summary(mod1)
summary(mod2)

# saprotrophs decreased with C:N, symbiotrophs increased with C:N, and NS for combined guild (p-0.06, +)

# step model for relative abundance of saprotrophs
all.parms <- lm(saprotroph ~ organic.soil.C.N + mineral.soil.C.N + GWC + white.spruce + paper.birch + dwarf.birch + moss + vaccinium + fern + kinnikinnick + 
                  heather + eriophorum + fireweed + horsetail + alder + hogsweed + willow + redberry + azalea + slope + ln(Distance) + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(formula = saprotroph ~ organic.soil.C.N + mineral.soil.C.N + 
                        paper.birch + vaccinium + fern + hogsweed + azalea + slope + 
                        ln(Distance) + carcass.deposition + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# saprotroph relative abundance increased with carcass deposition, distance, and decreased with organic soil C:N


mod <- lm(saprotroph ~ GWC + paper.birch + moss + fern + horsetail + alder + hogsweed, data=exploration.guild.RA)
summary(mod)

mod <- betareg(saprotroph_trans ~ C.N + GWC + paper.birch + moss + fern + horsetail + alder + hogsweed, data=exploration.guild.RA)
summary(mod)

# saprotrophs decreased with C:N, increased with horsetail and hogsweed

# step model for relative abundance of symbiotrophs
all.parms <- lm(symbiotroph ~ organic.soil.C.N + mineral.soil.C.N + GWC + white.spruce + paper.birch + dwarf.birch + moss + vaccinium + fern + kinnikinnick + 
                       heather + eriophorum + fireweed + horsetail + alder + hogsweed + willow + redberry + azalea + slope + ln(Distance) + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(formula = symbiotroph ~ paper.birch + moss + fern + horsetail + hogsweed + willow + azalea + ln(Distance) + carcass.deposition + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# symbiotrophs decreased with distance and carcass deposition

mod <- lm(symbiotroph ~ C.N + fern + horsetail + hogsweed + willow + azalea + Distance, data=exploration.guild.RA)
summary(mod)

mod <- betareg(symbiotroph_trans ~ C.N + fern + horsetail + hogsweed + willow + azalea + Distance, data=exploration.guild.RA)
summary(mod)

# symbiotrophs increased with C:N, also with horsetail and azalea, decreased with distance 

#########################################
### OVERALL BY EXPLORATION TYPE ###
#########################################

# examine data by stream
ggplot(exploration.guild.RA,aes(x=short.distance, y=Stream)) + geom_boxplot(alpha=0.2) + stat_summary(fun=mean, geom="point")
ggplot(exploration.guild.RA,aes(x=medium.distance, y=Stream)) + geom_boxplot(alpha=0.2) + stat_summary(fun=mean, geom="point")
ggplot(exploration.guild.RA,aes(x=long.distance, y=Stream)) + geom_boxplot(alpha=0.2) + stat_summary(fun=mean, geom="point")

# after examining relative abundance of guilds by stream, maybe include stream for long-distance

# non parametric kendall-theill regresson
mod1 <- mblm(short.distance ~ C.N, data=exploration.guild.RA)
mod2 <- mblm(medium.distance ~ C.N, data=exploration.guild.RA)
mod3 <- mblm(long.distance ~ C.N, data=exploration.guild.RA)
summary(mod1)
summary(mod2)
summary(mod3)

# short distance and medium-distance increased with C:N, and long distance decreased with C:N

# LM and GLM
mod1 <- lm(short.distance ~ C.N, data=exploration.guild.RA)
mod2 <- lm(medium.distance ~ C.N, data=exploration.guild.RA)
mod3 <- lm(long.distance ~ C.N, data=exploration.guild.RA)
summary(mod1)
summary(mod2)
summary(mod3)

# long distance decreased with C:N

# SHORT DISTANCE #

# step model for relative abundance of short distance types
all.parms <- lm(short.distance ~ organic.soil.C.N + mineral.soil.C.N + GWC + white.spruce + paper.birch + dwarf.birch + moss + vaccinium + fern + kinnikinnick + 
                  heather + eriophorum + fireweed + horsetail + alder + hogsweed + willow + redberry + azalea + slope + ln(Distance) + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(formula = short.distance ~ moss + fern + fireweed + horsetail + azalea + ln(Distance) + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# short distance decreased with distance

mod <- lm(short.distance ~ fern + fireweed + horsetail + azalea + Distance, data=exploration.guild.RA)
summary(mod)

mod <- betareg(short.distance ~ C.N + fern + fireweed + horsetail + azalea + Distance, data=exploration.guild.RA)
summary(mod)

# short.distance increased with azalea, and decreased with fireweed and distance

all.parms <- lm(medium.distance.fringe ~ organic.soil.C.N + mineral.soil.C.N + GWC + white.spruce + paper.birch + dwarf.birch + moss + vaccinium + fern + kinnikinnick + 
                  heather + eriophorum + fireweed + horsetail + alder + hogsweed + willow + redberry + azalea + slope + ln(Distance) + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(formula = medium.distance.fringe ~ organic.soil.C.N + mineral.soil.C.N + azalea + slope + carcass.deposition + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# medium distance fringe increased with slope and organic soil C:N, and decreased with mineral soil C:N
# without plants, decreased with carcass deposition and mineral soil C:N, and increased with slope and organic soil C:N

mod <- lm(medium.distance ~ C.N + paper.birch + fireweed + azalea + slope, data=exploration.guild.RA)
summary(mod)

mod <- betareg(medium.distance_trans ~ C.N + paper.birch + fireweed + azalea + slope, data=exploration.guild.RA)
summary(mod)

# medium distance increased with C:N, also increased with paper birch, slope, and decomposing carcass!

# MEDIUM DISTANCE SMOOTH # 

all.parms <- lm(medium.distance.smooth ~ organic.soil.C.N + mineral.soil.C.N + GWC + white.spruce + paper.birch + dwarf.birch + moss + vaccinium + fern + kinnikinnick + 
                  heather + eriophorum + fireweed + horsetail + alder + hogsweed + willow + redberry + azalea + slope + ln(Distance) + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(formula = medium.distance.smooth ~ paper.birch + dwarf.birch + fireweed + horsetail + azalea + ln(Distance) + decomposing.carcass + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# medium distance smooth increase with decomposing carcass and decreased with distance

# LONG DISTANCE # 

all.parms <- lm(long.distance ~ organic.soil.C.N + mineral.soil.C.N + GWC + white.spruce + paper.birch + dwarf.birch + moss + vaccinium + fern + kinnikinnick + 
                  heather + eriophorum + fireweed + horsetail + alder + hogsweed + willow + redberry + azalea + slope + ln(Distance) + 
                  carcass.deposition + decomposing.carcass, data=exploration.guild.RA)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

step.model.lmer <- lmer(formula = long.distance ~ paper.birch + moss + eriophorum + fireweed + ln(Distance) + decomposing.carcass + (1|Stream:Transect), data = exploration.guild.RA)
summary(step.model.lmer)

# long distance decreased with distance and increased with decomposing carcass
# without plants, long distance increased with decomposing carcass

mod <- lm(long.distance ~ white.spruce + paper.birch + moss + fern + decomposing.carcass, data=exploration.guild.RA)
summary(mod)

mod <- betareg(long.distance_trans ~ C.N, data=exploration.guild.RA)
summary(mod)
# this betareg doesn't work

# long distance increased with decomposing carcass, and decreased with paper birch presence


#########################################
### HANSEN OVERALL ####
#########################################

hansen$ph <- as.numeric(hansen$ph)

all.parms <- mod1 <- lm(saprotroph ~ C.N:Side + GWC + ph + slope + Side:Distance + carcass.deposition + decomposing.carcass, data=hansen)
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)
mod3 <- lm(saprotroph ~ C.N + GWC + slope + Side:Distance + carcass.deposition + decomposing.carcass, data=hansen)

mod1 <- lm(short.distance ~ C.N:Side, data=hansen)
mod2 <- lm(medium.distance ~ C.N:Side, data=hansen)
mod3 <- lm(long.distance ~ C.N + Side, data=hansen)
summary(mod1)
summary(mod2)
summary(mod3)

summary(mod1)
summary(mod2)
summary(mod3)
# no differences

#########################################
### HANSEN BY SIDE ####
#########################################

#########################################
### HANSEN LEFT ####
#########################################

ra.hansen <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Hansen"),]

mod1 <- mblm(saprotroph ~ C.N, data=ra.hansen)
mod2 <- lm(symbiotroph ~ C.N:Side, data=ra.hansen)
summary(mod1)
summary(mod2)

ra.hansen.left <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Hansen" & exploration.guild.RA$Side=="SL"),]
ra.hansen.left.close <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Hansen" & exploration.guild.RA$Side=="SL" & exploration.guild.RA$Distance < 9),]

mod1 <- mblm(short.distance ~ C.N, data=ra.hansen.left)
mod2 <- mblm(medium.distance ~ C.N, data=ra.hansen.left)
mod3 <- mblm(long.distance ~ C.N, data=ra.hansen.left.close)
summary(mod1)
summary(mod2)
summary(mod3)
# medium distance and long distance increased with C:N at SL

mod1 <- mblm(saprotroph ~ C.N, data=ra.hansen.left)
mod2 <- mblm(symbiotroph ~ C.N, data=ra.hansen.left)
mod3 <- mblm(saprotroph.symbiotroph ~ C.N, data=ra.hansen.left)
summary(mod1)
summary(mod2)
summary(mod3)
# no differences

#########################################
### HANSEN RIGHT ###
#########################################

ra.hansen.right <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Hansen" & exploration.guild.RA$Side=="SR"),]
ra.hansen.right.close <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Hansen" & exploration.guild.RA$Side=="SR" & exploration.guild.RA$Distance < 21),]

mod1 <- mblm(short.distance ~ C.N, data=ra.hansen.right.close)
mod2 <- mblm(medium.distance ~ C.N, data=ra.hansen.right.close)
mod3 <- mblm(long.distance ~ C.N, data=ra.hansen.right.close)
summary(mod1)
summary(mod2)
summary(mod3)
# medium distance increased with C:N at SR, long distance maybe decreased (p = 0.09)s

mod1 <- mblm(saprotroph ~ C.N, data=ra.hansen.right)
mod2 <- mblm(symbiotroph ~ C.N, data=ra.hansen.right)
mod3 <- mblm(saprotroph.symbiotroph ~ C.N, data=ra.hansen.right)
summary(mod1)
summary(mod2)
summary(mod3)
# symbiotrophs increased with C:N! 

#########################################
### HAPPY ###
#########################################

happy <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Happy"),]

mod1 <- mblm(short.distance ~ C.N, data=happy)
mod2 <- mblm(medium.distance ~ C.N, data=happy)
mod3 <- mblm(long.distance ~ C.N, data=happy)
summary(mod1)
summary(mod2)
summary(mod3)

# short distance increased with C.N, medium and long distance decreased with C:N

mod1 <- mblm(saprotroph ~ C.N, data=happy)
mod2 <- mblm(symbiotroph ~ C.N, data=happy)
mod3 <- mblm(saprotroph.symbiotroph ~ C.N, data=happy)
summary(mod1)
summary(mod2)
summary(mod3)

# saprotrophs decreased with C:N, symbiotrophs increased with C:N, combined guilds increased with C:N

#########################################
### YAKO ###
#########################################

yako <- exploration.guild.RA[which(exploration.guild.RA$Stream=="Yako"),]

mod1 <- mblm(short.distance ~ C.N, data=yako)
mod2 <- mblm(medium.distance ~ C.N, data=yako)
mod3 <- mblm(long.distance ~ C.N, data=yako)
summary(mod1)
summary(mod2)
summary(mod3)

# short distance increased with C.N, long distance decreased with C:N

mod1 <- mblm(saprotroph ~ C.N, data=yako)
mod2 <- mblm(symbiotroph ~ C.N, data=yako)
mod3 <- mblm(saprotroph.symbiotroph ~ C.N, data=yako)
summary(mod1)
summary(mod2)
summary(mod3)

# saprotrophs decreased with C:N, symbiotrophs increased with C:N

#######################################################################
# EXPLORATORY GRAPHS # 
#######################################################################

plot(ra.hansen.left$Distance, ra.hansen.left$long.distance)
plot(ra.hansen.left$Distance, ra.hansen.left$short.distance)
plot(ra.hansen.left$Distance, ra.hansen.left$medium.distance)

plot(ra.hansen.right$Distance, ra.hansen.right$long.distance)
plot(ra.hansen.right$Distance, ra.hansen.right$short.distance)
plot(ra.hansen.right$Distance, ra.hansen.right$medium.distance)

plot(happy$Distance, happy$long.distance)
plot(happy$Distance, happy$short.distance)
plot(happy$Distance, happy$medium.distance)

plot(yako$Distance, yako$long.distance)
plot(yako$Distance, yako$short.distance)
plot(yako$Distance, yako$medium.distance)

#######################################################################
# more exploratory figures # 
#######################################################################

hansen <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen"),]

# plot long distance 
left <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SL" & exploration.guild.RA$Distance < 11),]
right <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SR" & exploration.guild.RA$Distance < 11),]

summarySE(hansen, measurevar ="long.distance", groupvars=c("Side"))

# for right bank, group value from 5 into 6
right$Distance[which(right$Distance==5)] <- 6
tgc.right <- summarySE(right, measurevar="long.distance", groupvars=c("Distance"))
tgc.right

tgc.left <- summarySE(left, measurevar="long.distance", groupvars=c("Distance"))
tgc.left

tgc.left[3,5] <- 0.0014

long.left <- ggplot(tgc.left, aes(Distance, long.distance)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Salmon-enhanced")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=long.distance-se, ymax=long.distance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14),
        legend.position="none", plot.title=element_text(hjust=0.5)) + scale_x_reverse(breaks=c(1,3,6,10,20)) + ylim(0,0.0045) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) 

long.right <- ggplot(tgc.right, aes(Distance, long.distance)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() +ggtitle("Salmon-depleted")+
  theme_classic() + ylab("Relative Abundance") + geom_errorbar(aes(ymin=long.distance-se, ymax=long.distance+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_blank(), axis.text.y = element_blank(), axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14),
        legend.position="none", plot.title=element_text(hjust=0.5))  + 
  scale_x_continuous(breaks=c(1,3,6,10,20), expand = c(0.05, 0.05)) +  ylim(0,0.0045) + scale_y_continuous(expand = expansion(mult = c(0, 0.05)))

p <- grid.arrange(long.left, long.right, nrow=1, ncol=2)
ggsave(plot=p, width=10, height=6, dpi=300, filename = "long by bank.jpg")

left <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SL"),]
right <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SR"),]
tgc <- summarySE(left, measurevar="long.richness", groupvars=c("Distance"))
tgc

ggplot(tgc, aes(Distance, long.richness)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Long-distance types")+
  theme_classic() + ylab("Species richness") + geom_errorbar(aes(ymin=long.richness-se, ymax=long.richness+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14),
        legend.position="none") 

left <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SL"),]
right <- exploration.guild.RA[which(exploration.guild.RA$Stream == "Hansen" & exploration.guild.RA$Side=="SR"),]
tgc <- summarySE(right, measurevar="long.diversity", groupvars=c("Distance"))
tgc

ggplot(tgc, aes(Distance, long.diversity)) + geom_bar(position = "dodge", stat = "summary", fun = "mean") + theme_classic() + ggtitle("Long-distance types")+
  theme_classic() + ylab("Species richness") + geom_errorbar(aes(ymin=long.diversity-se, ymax=long.diversity+se), width=.2,position=position_dodge(.9)) + 
  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14),
        legend.position="none") 

#######################################################################
# final figures # 
#######################################################################

saprotroph <- ggplot(exploration.guild.RA, aes(C.N, saprotroph)) + geom_point() + stat_smooth(method="glm", method.args = list(family = "quasipoisson"), 
         formula = y ~ ns(x, 2)) + theme_classic() + ylab("Saprotroph relative abundance")+  theme(axis.title.y = element_text(size = 16), 
         axis.text.y = element_text(size = 14), axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none")

saprotroph.richness <- ggplot(exploration.guild.RA, aes(C.N, saprotroph.richness)) + geom_point() + stat_smooth(method="glm", method.args = list(family = "quasipoisson"), 
        formula = y ~ ns(x, 2)) + theme_classic() + ylab("Saprotroph species richness")+  theme(axis.title.y = element_text(size = 16), 
       axis.text.y = element_text(size = 14), axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none")

symbiotroph <- ggplot(exploration.guild.RA, aes(C.N, symbiotroph)) + geom_point() + geom_smooth(method='lm', formula = y~x) + theme_classic() + theme_classic() + 
  ylab("Symbiotroph relative abundance")+  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), 
  axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none")

symbiotroph.richness <- ggplot(exploration.guild.RA, aes(C.N, symbiotroph.richness)) + geom_point() + geom_smooth(method='lm', formula = y~x) + theme_classic() + theme_classic() + 
  ylab("Symbiotroph species richness")+  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), 
  axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none")

symbiotroph.diversity <- ggplot(exploration.guild.RA, aes(C.N, symbiotroph.diversity)) + geom_point() + geom_smooth(method='lm', formula = y~x) + theme_classic() + theme_classic() + 
  ylab("Symbiotroph species diversity")+  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), 
  axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none")

long.distance.all <- ggplot(exploration.guild.RA, aes(C.N, long.distance)) + geom_point() + stat_smooth(method="glm", 
                    method.args = list(family = "quasipoisson"), formula = y ~ ns(x, 1)) + ylim(0,0.06) + theme_classic() + 
  ylab("Relative Abundance of long distance types")+  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), 
  axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none") 

long.richness <- ggplot(exploration.guild.RA, aes(C.N, long.richness)) + geom_point() + stat_smooth(method="glm", 
   method.args = list(family = "quasipoisson"), formula = y ~ ns(x, 1)) + ylim(0,0.06) + theme_classic() + 
  ylab("Relative Abundance of long distance types")+  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), 
  axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none") 

short.distance.all <- ggplot(exploration.guild.RA, aes(C.N, short.distance)) + geom_point() + geom_smooth(method='lm', formula = y~x) + theme_classic() + theme_classic() + 
  ylab("Relative Abundance of short distance types")+  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), 
                                                    axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none")

short.distance.diversity <- ggplot(exploration.guild.RA, aes(C.N, short.diversity)) + geom_point() + geom_smooth(method='lm', formula = y~x) + theme_classic() + theme_classic() + 
  ylab("Species diversity of short distance types")+  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), 
  axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none")

medium.distance <- ggplot(exploration.guild.RA, aes(C.N, medium.distance)) + geom_point() + geom_smooth(method='lm', formula = y~x) + theme_classic() + theme_classic() + 
  ylab("Medium distance fungal type \n relative abundance")+  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), 
 axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none")

long.distance.hansen.left <- ggplot(ra.hansen.left.close, aes(C.N, long.distance)) + geom_point() + geom_smooth(method='lm', formula = y~x) + theme_classic() + theme_classic() + 
  ylab("Relative Abundance of long distance types")+  theme(axis.title.y = element_text(size = 16), axis.text.y = element_text(size = 14), 
   axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 14), legend.position="none")

p <- grid.arrange(saprotroph, symbiotroph, medium.distance, nrow=1, ncol=3)
setwd("~/Research/Mycorrhizae/Fieldwork/Genetics/NGS sequencing")
ggsave(plot=p, width=10, height=4, dpi=300, filename = "RA of guilds.jpg")

